void calc_orb(P *p, TRHO *rho, int *J, int nn, bool *failflag ) {
	
			if (8<nn) { 
				printf("ERROR: nn=%i, max is 8\n", nn);
				exit(0);
			}
	
			REAL pos[8][3], vel[8][3], dt[8], dt_max[8];
			REAL acc[8][3], t[8];

			bool resultF = true;
			int nsteps = ns*np;
			REAL ins2PI = 1.0/ns*2*M_PI;
			REAL eps2 = SQR(eps);
			
			
			for (int j=0; j<nn; j++)
			memset(&rho[J[j]*dg.n], 0, dg.n*sizeof(*rho));
			
			
			for (int j=0; j<nn; j++) {
				
			memcpy(pos[j], p[J[j]].Pos, 3*sizeof(REAL)); 
			memcpy(vel[j], p[J[j]].Vel, 3*sizeof(REAL));
			
			REAL r2 = SQR(pos[j][0])+SQR(pos[j][1]);
			REAL r = sqrt(r2);
			REAL z = fabs(pos[j][2]);
			
			// serial
			
			AG::A a = ag.get_acc_pot(r,z);
			
			
			
			REAL a_mod = sqrt( SQR(a.r)+SQR(a.z) );
			
			if (r!=0) {
				REAL a_r_over_r(a.r/r);
				acc[j][0]= pos[j][0]*a_r_over_r;
				acc[j][1]= pos[j][1]*a_r_over_r;
			} 
			
			if (z!=0)
				acc[j][2]= pos[j][2]*a.z/z;
	
			
			REAL posm = sqrt( SQR(pos[j][0]) + SQR(pos[j][1]) + SQR(pos[j][2]) );
			REAL velm = sqrt( SQR(vel[j][0]) + SQR(vel[j][1]) + SQR(vel[j][2]) );
			REAL accm = sqrt( SQR(acc[j][0]) + SQR(acc[j][1]) + SQR(acc[j][2]) );
			
			
			REAL tp = 2*M_PI*posm/sqrt(accm*posm);
			dt[j] = tp/ns;

			
			t[j] = 0;
			
			
			dt_max[j] = 0.1*2*M_PI*r*sqrt(-a.pot);
			
			// 0. predict velocity at dt/2
			vel[j][0] += acc[j][0]*.5*dt[j];
			vel[j][1] += acc[j][1]*.5*dt[j];
			vel[j][2] += acc[j][2]*.5*dt[j];
		
			failflag[j] = false;
				
			}
			
			
 			for (int i=0; i<nsteps; i++) {
			for (int j=0; j<nn; j++) {
				// 1. predict new position 
				pos[j][0] += vel[j][0]*dt[j];
				pos[j][1] += vel[j][1]*dt[j];
				pos[j][2] += vel[j][2]*dt[j];
				
				// 2. get accelerations
				REAL r2 = SQR(pos[j][0])+SQR(pos[j][1]);
				REAL r = sqrt(r2);
				REAL z = fabs(pos[j][2]);
				
				REAL R2 = r2 + SQR(pos[j][2]);
				REAL R = sqrt(R2);
				
				AG::A  a = ag.get_acc(r,z);
				
					
				REAL a_mod = sqrtf(SQR(a.r)+SQR(a.z));
				REAL dt_new = (ins2PI*R)/sqrtf(a_mod*R);
				
				
				
				// 2. get accelerations
				
				#ifdef DG_DTMAX_CELL
				REAL vr = fabs(vel[j][0]*pos[j][0] + vel[j][1]*pos[j][1])/r;
				REAL vz = fabs(vel[j][2]);
				
				int ir = dg.r2ir(r);
				int iz = dg.z2iz(z);
				
				REAL dr = dg.ir2r(ir+1)-dg.ir2r(ir);
				REAL dz = dg.iz2z(iz+1)-dg.iz2z(iz);
				
				REAL dtrmax = dr/vr;
				REAL dtzmax = dz/vz;
				
				if (dtrmax<dt_new) dt_new = dtrmax; 
				if (dtzmax<dt_new) dt_new = dtzmax;
				#endif
				
				
				if (a.r==0 && a.z==0) failflag[j]=true;
				
				if (dt_max[j]<dt_new) dt_new = dt_max[j];
				
				REAL _dt = 0.5*(dt[j] + dt_new);
				t[j] +=_dt;
				dt[j] = dt_new;
				
				//if (!failflag[j]) 
				#ifdef DG_ADD_STD
				dg.add(&rho[J[j]*dg.n], r, z, _dt);
				#endif
				

				
				
				
				if (r!=0) {
					REAL a_r_over_r(a.r/r);
					acc[j][0]= pos[j][0]*a_r_over_r;
					acc[j][1]= pos[j][1]*a_r_over_r;
				} 
			
				if (z!=0)
					acc[j][2]= pos[j][2]*(a.z/z);
		
				/*
				if (j==0) {
					//a.print();
					printf("%4i % .7e, % .7e\n", i, r, z);
					//printf("%i %g, %g\n", i, r, z);
				//if (r[0]!=r[0]) exit(0);
				}
				*/
			
				
				//printf("%i %g, %g\n", i, r,z);
				
				
				// 3. predic new velocities				
				vel[j][0] += acc[j][0]*_dt;
				vel[j][1] += acc[j][1]*_dt;
				vel[j][2] += acc[j][2]*_dt;
				
				//printf("vel = %g %g %g\n\n", vel[j][0], vel[j][1], vel[j][2]);
				/*
				if (i==1000) {
					printf("cala_orb()::exit(0);\n");
					exit(0);
				}
				*/
			}
			}
			
			for (int j=0; j<nn; j++) {
				pos[j][0] += vel[j][0]*0.5*dt[j];
				pos[j][1] += vel[j][1]*0.5*dt[j];
				pos[j][2] += vel[j][2]*0.5*dt[j];
			}
			
			
			for (int j=0; j<nn; j++)
				if (t[j]>0) {
					REAL m_over_t(p[J[j]].Mass/t[j]);
					REAL *_rho = &rho[J[j]*dg.n];
					for(int i=0; i<dg.n; i++)  
					if (_rho[i]) _rho[i]*=m_over_t;
				}
				
			
			//return resultF;
} 
